{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8de3639f-7e75-4d04-81e5-6827dc81b72a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np\n",
    "from toolkits import gauss_line_int, gauss_tri_int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "cfb90fd5-ba63-4daf-8d99-f688b3d4a1ef",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class Poisson2D:\n",
    "    \"\"\"docstring for Poisson\n",
    "    \n",
    "    3-Node triangular element FEM/EFEM model for 2D Poisson/Laplace equation. \n",
    "\n",
    "    \"\"\"\n",
    "    def __init__(self, geometry, DirMarker, NeuMarker, holeMarker, BCfunc):\n",
    "        self.Geometry = geometry\n",
    "\n",
    "        self.Node = None\n",
    "        self.Edge = dict()\n",
    "        self.Element = None\n",
    "\n",
    "        self.DirMarker = DirMarker # list of markers for Dirichlet's boundary\n",
    "        self.Dirichlet = dict()    # (i,j): marker\n",
    "        # self.DirEdge = dict()\n",
    "        self.NeuMarker = NeuMarker # list of markers for Neumann's boundary\n",
    "        self.Neumann = dict()      # (i,j): marker\n",
    "        # self.NeuEdge = dict()\n",
    "        self.NeuCluster = None\n",
    "        self.holeMarker = holeMarker  # list of marker lists for each hole \n",
    "        self.hole = dict()      # (i,j): holeid \n",
    "        self.holeEdge = dict()  # holeid : [(i,j),(m,n),...]\n",
    "\n",
    "        self.crossEdge = dict()\n",
    "        self.crossElement = dict()\n",
    "\n",
    "        self.BCfunc = BCfunc\n",
    "\n",
    "        self.numHole = len(self.holeMarker)\n",
    "        self.numElement = None\n",
    "        self.numNode = None\n",
    "        self.numDOF = None\n",
    "\n",
    "        self.area = None\n",
    "        # EFEM matrix vector\n",
    "        self.FU = None\n",
    "        self.dU = None\n",
    "        self.F = None\n",
    "        self.d = None\n",
    "        self.PhiFe = None\n",
    "        self.s = None\n",
    "        self.qh = None\n",
    "        # FEM matrix vector\n",
    "\n",
    "        self._mesh = None\n",
    "        self._edof = None\n",
    "        self._BdElements = None\n",
    "      \n",
    "\n",
    "    def mesh(self, showMesh=False):\n",
    "        # generate mesh using GMSH.\n",
    "        self._mesh = cfm.GmshMesh(self.Geometry, return_boundary_elements=True)\n",
    "        self._mesh.elType = 2\n",
    "        self._mesh.dofsPerNode = 1\n",
    "        ##########################################\n",
    "        #                                        #\n",
    "        #          Basic Data Structure          #\n",
    "        #                                        #\n",
    "        # Node: list of coordinates (x, y)       #\n",
    "        # Element: list of nodes (i,j,k) CCW.    #\n",
    "        # Edge: dict of edges (i,j):n pair       #\n",
    "        #                                        #\n",
    "        ##########################################\n",
    "        self.Node, self._edof, _, _, _, self._BdElements = self._mesh.create()\n",
    "        self.Element = self._edof - 1\n",
    "        for n, e in enumerate(self.Element):\n",
    "            self.Edge[(e[0], e[1])] = n\n",
    "            self.Edge[(e[1], e[2])] = n\n",
    "            self.Edge[(e[2], e[0])] = n\n",
    "        self.numElement = len(self.Element)\n",
    "        self.numNode = len(self.Node)\n",
    "        self.numDOF = self.numNode + self.numHole\n",
    "        ###############################################\n",
    "        #                                             #\n",
    "        #                 Marker Bank                 #\n",
    "        #              (for applying BC)              #\n",
    "        #                                             #\n",
    "        # DirEdge: dict of marker, marker:(i,j)       #\n",
    "        # NeuEdge: dict of marker, marker:(i,j)       #\n",
    "        # holeEdge: dict of holeid, holeid:[(i,j),..] #\n",
    "        #                                             #\n",
    "        ###############################################\n",
    "        # for marker in self.DirMarker:\n",
    "        #     self.DirEdge[marker] = np.array([e['node-number-list'] for e in self._BdElements[marker]])-1\n",
    "        \n",
    "        # for marker in self.NeuMarker:\n",
    "        #     self.NeuEdge[marker] = np.array([e['node-number-list'] for e in self._BdElements[marker]])-1\n",
    "\n",
    "        for i, hole in enumerate(self.holeMarker):\n",
    "            edges = []\n",
    "            for marker in hole:\n",
    "                for e in self._BdElements[marker]:\n",
    "                    edge = tuple(np.array(e['node-number-list'])-1)\n",
    "                    # Force to CCW.\n",
    "                    if edge not in self.Edge:\n",
    "                        edge = edge[::-1]\n",
    "                    edges.append(edge)\n",
    "            self.holeEdge[i] = np.array(edges)\n",
    "\n",
    "        ##############################################\n",
    "        #                                            #\n",
    "        #                Boundary sets               #\n",
    "        #          (for applying constraint)         #\n",
    "        #                                            #\n",
    "        # Dirichlet: Dirichlet edge set.(i,j):marker #\n",
    "        # Neumann: Neumann edge set. (i,j):marker    #\n",
    "        # hole: hole edge dict. (i,j):holeid         #  \n",
    "        #                                            #\n",
    "        ##############################################\n",
    "        for marker in self.DirMarker:\n",
    "            for e in self._BdElements[marker]:\n",
    "                edge = tuple(np.array(e['node-number-list'])-1)\n",
    "                # Force to CCW.\n",
    "                if edge not in self.Edge:\n",
    "                    edge = edge[::-1]\n",
    "                self.Dirichlet[edge] = marker\n",
    "\n",
    "        for marker in self.NeuMarker:\n",
    "            for e in self._BdElements[marker]:\n",
    "                edge = tuple(np.array(e['node-number-list'])-1)\n",
    "                # Force to CCW.\n",
    "                if edge not in self.Edge:\n",
    "                    edge = edge[::-1]\n",
    "                self.Neumann[edge] = marker                \n",
    "\n",
    "        self.NeuCluster = self._cluster_edge(self.Neumann)\n",
    "\n",
    "        # self._check_ccw()\n",
    "\n",
    "        if len(self.holeMarker) != 0: \n",
    "            for n, hms in enumerate(self.holeMarker):\n",
    "                for m in hms:\n",
    "                    for e in self._BdElements[m]:\n",
    "                        edge = tuple(np.array(e['node-number-list'])-1)\n",
    "                        # Force to CCW.\n",
    "                        if edge not in self.Edge:\n",
    "                            edge = edge[::-1]\n",
    "                        self.hole[edge] = n\n",
    "            self._connect_domain()\n",
    "        # Display mesh\n",
    "        if showMesh:\n",
    "            cfv.draw_mesh(coords=self.Node, edof=self._edof, dofs_per_node=self._mesh.dofs_per_node,\n",
    "                          el_type=self._mesh.el_type,filled=True,title=\"Mesh\")\n",
    "            cfv.showAndWait()\n",
    "\n",
    "    def _cluster_edge(self, edges):\n",
    "        from collections import deque\n",
    "        from copy import deepcopy\n",
    "\n",
    "        edges = deepcopy(edges)\n",
    "        clusters = []\n",
    "\n",
    "        e, _ = edges.popitem()\n",
    "        clusters.append(deque([e]))\n",
    "        \n",
    "        while len(edges) != 0:\n",
    "            quit = False\n",
    "            for i, j in edges:\n",
    "                for c in clusters:\n",
    "                    if i == c[-1][-1]:\n",
    "                        c.append((i, j))\n",
    "                        edges.pop((i, j), None)\n",
    "                        quit = True\n",
    "                        break\n",
    "                    if j == c[0][0]:\n",
    "                        c.appendleft((i, j))\n",
    "                        edges.pop((i, j), None)\n",
    "                        quit = True\n",
    "                        break\n",
    "                if quit:\n",
    "                    break\n",
    "            else:\n",
    "                if len(edges) != 0:\n",
    "                    e, _ = edges.popitem()\n",
    "                    clusters.append(deque([e]))\n",
    "                    \n",
    "        # convert to list\n",
    "        clusters = [list(q) for q in clusters]\n",
    "        return clusters \n",
    "\n",
    "\n",
    "    def _check_ccw(self):\n",
    "        import matplotlib.pyplot as plt\n",
    "        plt.figure(figsize=(20,10),dpi=300)\n",
    "        plt.scatter(self.Node[:,0], self.Node[:,1])\n",
    "        for i,j in self.Neumann:\n",
    "            x1,y1 = self.Node[i]\n",
    "            x2,y2 = self.Node[j]\n",
    "            plt.arrow(x1,y2,x2-x1,y2-y1)\n",
    "        plt.savefig('ccw.jpg')\n",
    "        plt.show()\n",
    "\n",
    "\n",
    "    def _connect_domain(self):\n",
    "        #################################\n",
    "        #                               #\n",
    "        #  Crossing elements and edges  #\n",
    "        #                               #\n",
    "        # crossEdge: (i,j):(n,eta)      #\n",
    "        #                               #\n",
    "        #################################\n",
    "\n",
    "        direction = np.array([0, -1])\n",
    "\n",
    "        for n in self.holeEdge:\n",
    "            edge = tuple(self.holeEdge[n][0])\n",
    "            self.crossEdge[edge] = np.array([[n+self.numNode], [-1]])\n",
    "\n",
    "            _connected = False\n",
    "\n",
    "            while not _connected:\n",
    "                # print('current:', edge, self.Node[list(edge)][0], self.Node[list(edge)][1])\n",
    "                ei = self.Edge[edge]\n",
    "                \n",
    "                i, j, k = self.Element[ei]\n",
    "                order = [i,j,k]  \n",
    "                \n",
    "                local_edges = set([(i,j), (j,k), (k,i)])\n",
    "                local_edges.remove(edge)\n",
    "                e1 = np.array(local_edges.pop())\n",
    "                e2 = np.array(local_edges.pop())\n",
    "                \n",
    "                e1n = self.Node[e1[1]] - self.Node[e1[0]]\n",
    "                e1n = e1n/np.linalg.norm(e1n)\n",
    "                \n",
    "                e2n = self.Node[e2[1]] - self.Node[e2[0]]\n",
    "                e2n = e2n/np.linalg.norm(e2n)\n",
    "                \n",
    "                d1 = np.abs(np.dot(direction, e1n))\n",
    "                d2 = np.abs(np.dot(direction, e2n))\n",
    "                ae = np.zeros((3,1))\n",
    "                if d1 < d2:\n",
    "                    L1 = np.linalg.norm(self.Node[edge[0]] - self.Node[edge[1]])\n",
    "                    L2 = np.linalg.norm(self.Node[e1[0]] - self.Node[e1[1]])\n",
    "                    ae[order.index(edge[0])]= 1/L1\n",
    "                    ae[order.index(e1[0])]= -1/L2\n",
    "                    edge = tuple(e1)\n",
    "                else:\n",
    "                    L1 = np.linalg.norm(self.Node[edge[0]] - self.Node[edge[1]])\n",
    "                    L2 = np.linalg.norm(self.Node[e2[0]] - self.Node[e2[1]])\n",
    "                    ae[order.index(edge[0])]= 1/L1\n",
    "                    ae[order.index(e2[0])]= -1/L2\n",
    "                    edge = tuple(e2)\n",
    "\n",
    "                ae = np.insert(ae, 0, n+self.numNode, axis=0)\n",
    "                if ei in self.crossElement:\n",
    "                    self.crossElement[ei] = np.c_[self.crossElement[ei], ae]\n",
    "                else:\n",
    "                    self.crossElement[ei] = ae\n",
    "                # print(e1,'\\t',e1n)\n",
    "                # print(e2,'\\t',e2n)\n",
    "                # print(ae)\n",
    "                # print(d1,d2)\n",
    "                # print(\"----------------------------------\")\n",
    "                # print('next:', edge, self.Node[list(edge)][0], self.Node[list(edge)][1])\n",
    "                # print(self.Neumann.keys())\n",
    "                if (edge in self.Dirichlet) or (edge in self.Neumann):\n",
    "                    if edge in self.hole:\n",
    "                        m = self.hole[edge]\n",
    "                        if m == n:\n",
    "                            # TODO: reject.\n",
    "                            pass\n",
    "                        else:\n",
    "                            if edge in self.crossEdge:\n",
    "                                self.crossEdge[edge] = np.c_[self.crossEdge[edge], np.array([[n+self.numNode], [-1]])]\n",
    "                            else:\n",
    "                                self.crossEdge[edge] = np.array([[n+self.numNode], [-1]])\n",
    "                    else:\n",
    "                        if edge in self.crossEdge:\n",
    "                            self.crossEdge[edge] = np.c_[self.crossEdge[edge], np.array([[n+self.numNode], [1]])]\n",
    "                        else:\n",
    "                            self.crossEdge[edge] = np.array([[n+self.numNode], [1]])\n",
    "\n",
    "                    _connected = True\n",
    "                else:\n",
    "                    edge = edge[::-1]\n",
    "                    # print(self.Node[list(edge)])\n",
    "\n",
    "    def _assembly_efem(self):\n",
    "        ######################################\n",
    "        #                                    #\n",
    "        #  unconstrained flexibility matrix  #\n",
    "        #                                    #\n",
    "        ######################################\n",
    "        self.FU = np.zeros((self.numDOF, self.numDOF))\n",
    "        self.PhiFe = np.zeros((self.numElement, 2, 3))\n",
    "        self.area = np.zeros(self.numElement)\n",
    "\n",
    "        for i, e in enumerate(self.Element):\n",
    "            x1, y1 = self.Node[e][0]\n",
    "            x2, y2 = self.Node[e][1]\n",
    "            x3, y3 = self.Node[e][2]\n",
    "            \n",
    "            self.area[i] = 1/2*np.linalg.det(np.array([[x1, y1, 1],\n",
    "                                                       [x2, y2, 1],\n",
    "                                                       [x3, y3, 1]]))\n",
    "            \n",
    "            self.PhiFe[i] = np.array([[x2-x3, x3-x1, x1-x2],\n",
    "                                      [y2-y3, y3-y1, y1-y2]])/(2*self.area[i])\n",
    "            \n",
    "\n",
    "            if i in self.crossElement:\n",
    "                L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "                L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "                left = np.linalg.inv([[(y1-y2)/L1, (x2-x1)/L1],\n",
    "                                      [(y2-y3)/L2, (x3-x2)/L2]])\n",
    "\n",
    "                PhiLe = np.c_[left, np.zeros((2,1))]\n",
    "                \n",
    "                # print(Eh[i])\n",
    "                PhiFe = np.c_[self.PhiFe[i], PhiLe @ self.crossElement[i][1::,:]]\n",
    "                Fe = self.area[i] * PhiFe.T @ PhiFe\n",
    "                \n",
    "                ### Todo: multiple holes ###\n",
    "                ey = np.r_[e, self.crossElement[i][0,:]].astype(int)\n",
    "                self.FU[np.ix_(ey,ey)] = self.FU[np.ix_(ey,ey)] + Fe\n",
    "            else:\n",
    "                Fe = self.area[i] * self.PhiFe[i].T @ self.PhiFe[i]\n",
    "                \n",
    "                self.FU[np.ix_(e,e)] = self.FU[np.ix_(e,e)] + Fe\n",
    "\n",
    "        print('size: {} x {}'.format(*self.FU.shape))\n",
    "        print('rank: {}'.format(np.linalg.matrix_rank(self.FU)))\n",
    "\n",
    "        ######################################\n",
    "        #                                    #\n",
    "        #   unconstrained kinematic vector   #\n",
    "        #                                    #\n",
    "        ######################################\n",
    "        self.dU = np.zeros(self.numDOF)\n",
    "        #========================#\n",
    "        #   Dirichlet boundary   #\n",
    "        #========================#\n",
    "        for edge in self.Dirichlet:\n",
    "            marker = self.Dirichlet[edge]\n",
    "            i, j = edge\n",
    "            L = np.linalg.norm(self.Node[i]-self.Node[j])\n",
    "            up = gauss_line_int(self.BCfunc[marker], *self.Node[i], *self.Node[j])\n",
    "            self.dU[i] = self.dU[i] - up/L\n",
    "            self.dU[j] = self.dU[j] + up/L\n",
    "\n",
    "        # for marker in self.DirEdge:\n",
    "        #     for i,j in self.DirEdge[marker]:\n",
    "        #         L = np.linalg.norm(self.Node[i]-self.Node[j])\n",
    "        #         up = gauss_line_int(self.BCfunc[marker], *self.Node[i], *self.Node[j])\n",
    "        #         self.dU[i] = self.dU[i] - up/L\n",
    "        #         self.dU[j] = self.dU[j] + up/L\n",
    "\n",
    "\n",
    "    def _constraint(self):\n",
    "        #========================================#\n",
    "        #    Neumann boundary (without source)   #\n",
    "        #========================================#\n",
    "        self.sN = np.zeros(self.numDOF)\n",
    "        self.A = np.eye(self.numDOF)\n",
    "\n",
    "        # todo : move to subfunction.\n",
    "        eff_dofs = []\n",
    "        red_dofs = []\n",
    "        for c in self.NeuCluster:\n",
    "            eff_dofs.append([c[0][0]])\n",
    "            _red = []\n",
    "            for e in c:\n",
    "                # to avoid error occured in closed-edge.\n",
    "                if e[1] != eff_dofs[-1][0]:\n",
    "                    _red.append(e[1])\n",
    "            red_dofs.append(_red)\n",
    "\n",
    "        for n, d in enumerate(eff_dofs):\n",
    "            self.A[np.ix_(red_dofs[n], d)] = 1\n",
    "\n",
    "        all_red_dofs = [y for x in red_dofs for y in x]\n",
    "\n",
    "        # Neumann terms\n",
    "        for c in self.NeuCluster:\n",
    "            print('first Neumann edge:', self.Node[c[0][0]], self.Node[c[0][1]])\n",
    "            rn = {}\n",
    "            fl_bar = 0\n",
    "            # exclude the last edge in closed edge.\n",
    "            if c[0][0] == c[-1][-1]:\n",
    "                cn = c[:-1]\n",
    "            else:\n",
    "                cn = c\n",
    "\n",
    "            for e in cn:\n",
    "                marker = self.Neumann[e]\n",
    "                i, j = e    \n",
    "                fl_bar = fl_bar + gauss_line_int(self.BCfunc[marker], *self.Node[i], *self.Node[j])\n",
    "                self.sN[j] = fl_bar\n",
    "\n",
    "                if e in self.crossEdge:\n",
    "                    print(e, self.Node[e[0]], self.Node[e[1]])\n",
    "                    dofs, etas = self.crossEdge[e]\n",
    "                    for dof, eta in zip(dofs, etas):\n",
    "                        rn[int(dof)] = eta\n",
    "\n",
    "                for dof in rn:\n",
    "                    self.A[j, dof] = rn[dof]\n",
    "\n",
    "        self.A = np.delete(self.A, all_red_dofs, axis=1)\n",
    "\n",
    "        # rk \n",
    "\n",
    "        FU = self.FU[:]\n",
    "        dU = self.dU[:]\n",
    "\n",
    "        for c in self.NeuCluster:\n",
    "            # 封闭的一定是孔吗？\n",
    "            if c[0][0] == c[-1][-1]:\n",
    "                n = self.hole[c[0]]\n",
    "                dof = n + self.numNode\n",
    "                r = 0\n",
    "                for e in c:\n",
    "                    marker = self.Neumann[e]\n",
    "                    r = r + gauss_line_int(self.BCfunc[marker], *self.Node[i], *self.Node[j])\n",
    "\n",
    "                FU[dof,:] = 0\n",
    "                FU[dof,dof] = 1\n",
    "                dU[dof] = r\n",
    "\n",
    "        # r1 = 0\n",
    "        FU[0,:] = 0\n",
    "        FU[0,0] = 1\n",
    "        dU[0] = 0\n",
    "\n",
    "        self.F = self.A.T @ FU @ self.A\n",
    "        self.d = self.A.T @ (dU - FU @ self.sN)\n",
    "\n",
    "\n",
    "    def _solve_efem(self):\n",
    "        s = np.linalg.inv(self.F) @ self.d\n",
    "        self.s = self.A @ s + self.sN\n",
    "\n",
    "    def _construct_flux(self):\n",
    "        self.qh = np.zeros((self.numElement, 2)) \n",
    "\n",
    "        for i, e in enumerate(self.Element):\n",
    "            \n",
    "            if i in self.crossElement:\n",
    "                x1,y1 = self.Node[e][0]\n",
    "                x2,y2 = self.Node[e][1]\n",
    "                x3,y3 = self.Node[e][2]\n",
    "                L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "                L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "                left = np.linalg.inv([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                                      [(y2-y3)/L2,(x3-x2)/L2]])\n",
    "\n",
    "                PhiLe = np.c_[left, np.zeros((2,1))]\n",
    "                \n",
    "                PhiFe = np.c_[self.PhiFe[i], PhiLe @ self.crossElement[i][1::,:]]\n",
    "\n",
    "                ey = np.r_[e, self.crossElement[i][0,:]].astype(int)\n",
    "\n",
    "                self.qh[i] = PhiFe @ self.s[ey]\n",
    "\n",
    "            else:\n",
    "                self.qh[i] = self.PhiFe[i] @ self.s[e]\n",
    "\n",
    "    def _display_efem(self):\n",
    "        cfv.figure()\n",
    "\n",
    "        # Draw the mesh.\n",
    "\n",
    "        cfv.drawElementValues(\n",
    "            ev=self.qh[:,1],\n",
    "            coords=self.Node,\n",
    "            edof=self._edof,\n",
    "            dofs_per_node=self._mesh.dofsPerNode,\n",
    "            el_type=self._mesh.elType,\n",
    "            title=\"EFEM flux field\",\n",
    "            draw_elements=True\n",
    "                )\n",
    "        cfv.showAndWait()\n",
    "\n",
    "\n",
    "    def efem(self):\n",
    "        if self._mesh == None:\n",
    "            print('Please use .meshing() method to create mesh at first!')\n",
    "        else:\n",
    "            self._assembly_efem()\n",
    "            self._constraint()\n",
    "            self._solve_efem()\n",
    "            self._construct_flux()\n",
    "            self._display_efem()\n",
    "\n",
    "    def fem(self):\n",
    "        if self._mesh == None:\n",
    "            print('Please use .meshing() method to create mesh at first!')\n",
    "        else:\n",
    "            pass\n",
    "\n",
    "    def cre(self):\n",
    "        pass\n",
    "\n",
    "    def diplay(self):\n",
    "        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7682be08-62ce-40f9-abc3-5480b5c12b87",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "size: 1311 x 1311\n",
      "rank: 1310\n",
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x0000015A970E0AC0>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def create_geometry(el_on_curve=5):\n",
    "    # Geometry definition\n",
    "    g = cfg.Geometry()\n",
    "\n",
    "    r1 = 0.15\n",
    "    r2 = r1/np.sqrt(2)\n",
    "\n",
    "    points = [[0,1],[0.4,1],[0.6,1],[0.8,1],[1.2,1],[1.4,1],[1.6,1],[2.0,1],  # 0-7\n",
    "              [0,0.7],[0.4,0.7],[0.6,0.7],[0.8,0.7],[1.2,0.7],[1.4,0.7],[1.6,0.7],[2.0,0.7], # 8-15\n",
    "              [0,0.5],[0.4,0.5],[0.6-r1,0.5],[0.6-r2,0.5+r2],[0.6,0.5+r1],[0.6+r2,0.5+r2],[0.6+r1,0.5],[0.8,0.5], # 16-23\n",
    "              [1.2,0.5],[1.4-r1,0.5],[1.4-r2,0.5+r2],[1.4,0.5+r1],[1.4+r2,0.5+r2],[1.4+r1,0.5],[1.6,0.5],[2,0.5], # 24-31\n",
    "              [0.6,0.5],[1.4,0.5], # 32-33 circle center\n",
    "              [0.6-r2,0.5-r2],[0.6,0.5-r1],[0.6+r2,0.5-r2],[1.4-r2,0.5-r2],[1.4,0.5-r1],[1.4+r2,0.5-r2], # 34-39\n",
    "              [0,0.3],[0.4,0.3],[0.6,0.3],[0.8,0.3],[1.2,0.3],[1.4,0.3],[1.6,0.3],[2,0.3], # 40-47\n",
    "              [0,0],[0.4,0],[0.6,0],[0.8,0],[1.2,0],[1.4,0],[1.6,0],[2,0]] # 48-55\n",
    "\n",
    "    splines = [[0,1],[1,2],[2,3],[3,4],[4,5],[5,6],[6,7], # 0-6 top\n",
    "               [0,8],[1,9],[2,10],[3,11],[4,12],[5,13],[6,14],[7,15], # 7-14\n",
    "               [8,9],[9,10],[10,11],[11,12],[12,13],[13,14],[14,15], # 15-21\n",
    "               [8,16],[9,17],[9,19],[10,20],[11,21],[11,23],[12,24],[12,26],[13,27],[14,28],[14,30],[15,31], # 22-33\n",
    "               [16,17],[17,18],[22,23],[23,24],[24,25],[29,30],[30,31], #34-40\n",
    "               [16,40],[17,41],[34,41],[35,42],[36,43],[23,43],[24,44],[37,44],[38,45],[39,46],[30,46],[31,47], #41-52\n",
    "               [40,41],[41,42],[42,43],[43,44],[44,45],[45,46],[46,47], # 53-59\n",
    "               [40,48],[41,49],[42,50],[43,51],[44,52],[45,53],[46,54],[47,55],# 60-67\n",
    "               [48,49],[49,50],[50,51],[51,52],[52,53],[53,54],[54,55],# 68-74\n",
    "              ]\n",
    "\n",
    "    longs = [0,3,6,15,18,21,34,37,40,53,56,59,68,71,74]\n",
    "\n",
    "    arcs = [[18,32,19],[19,32,20],[20,32,21],[21,32,22],\n",
    "            [22,32,36],[36,32,35],[35,32,34],[34,32,18], # 75-82\n",
    "            [25,33,26],[26,33,27],[27,33,28],[28,33,29],\n",
    "            [29,33,39],[39,33,38],[38,33,37],[37,33,25], # 83-90\n",
    "           ]\n",
    "\n",
    "    surfaces = [[0,8,15,7],[1,9,16,8],[2,10,17,9],[3,11,18,10],[4,12,19,11],[5,13,20,12],[6,14,21,13],\n",
    "             [15,23,34,22],[24,75,35,23],[16,25,76,24],[17,26,77,25],[26,27,36,78],\n",
    "             [18,28,37,27],\n",
    "             [29,83,38,28],[19,30,84,29],[20,31,85,30],[31,32,39,86],[21,33,40,32],\n",
    "             [34,42,53,41],[35,82,43,42],[81,44,54,43],[80,45,55,44],[36,46,45,79],\n",
    "             [37,47,56,46],\n",
    "             [38,90,48,47],[89,49,57,48],[88,50,58,49],[39,51,50,87],[40,52,59,51],\n",
    "             [53,61,68,60],[54,62,69,61],[55,63,70,62],[56,64,71,63],[57,65,72,64],[58,66,73,65],[59,67,74,66]]\n",
    "\n",
    "    for x, y in points:\n",
    "        g.point([x, y])\n",
    "\n",
    "    for i, s in enumerate(splines):\n",
    "        if i in longs:\n",
    "            g.spline(s, el_on_curve=2*el_on_curve)\n",
    "        else:\n",
    "            g.spline(s, el_on_curve=el_on_curve)\n",
    "\n",
    "    for c in arcs:\n",
    "        g.circle(c, el_on_curve=el_on_curve)\n",
    "\n",
    "    for f in surfaces:\n",
    "        g.struct_surf(f)\n",
    "\n",
    "\n",
    "    bcs = {100:[7,22,41,60],\n",
    "           200:[68,69,70,71,72,73,74],\n",
    "           300:[75,76,77,78],\n",
    "           101:[0,1,2,3,4,5,6],\n",
    "           201:[14,33,52,67],\n",
    "           301:[79,80,81,82],\n",
    "           401:[83,84,85,86,87,88,89,90],\n",
    "           }\n",
    "\n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)\n",
    "\n",
    "    return g\n",
    "\n",
    "def dir100(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir200(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir300(x,y):\n",
    "    return np.ones_like(x)\n",
    "\n",
    "def neu101(x,y):\n",
    "    return np.ones_like(x)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return np.ones_like(x)\n",
    "\n",
    "def neu301(x,y):\n",
    "    return np.ones_like(x)\n",
    "\n",
    "def neu401(x,y):\n",
    "    return np.ones_like(x)\n",
    "\n",
    "\n",
    "from fem import Poisson2D\n",
    "geometry = create_geometry(el_on_curve=5)\n",
    "Dirichlet = [100,200,300]\n",
    "Neumann = [101,201,301,401]\n",
    "hole = [[300,301],[401]]\n",
    "\n",
    "BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "          101:neu101, 201:neu201, 301:neu301, 401:neu401}\n",
    "\n",
    "\n",
    "model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc, direction=[[1,0],[-np.sqrt(2)/2,-np.sqrt(2)/2]])\n",
    "model.mesh()\n",
    "model.efem()\n",
    "# model.fem()\n",
    "\n",
    "# model.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7cad5b6e-5f7c-45dc-baf4-ad0989f012a9",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "size: 661 x 661\n",
      "rank: 660\n",
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x0000022FB17E6EB0>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def create_geometry(el_on_curve=5):\n",
    "    # Geometry definition\n",
    "    g = cfg.Geometry()\n",
    "\n",
    "    points = [[0,1],[0.4,1], [0.6,1], [1,1],\n",
    "              [0,0.6], [0.4,0.6], [0.6,0.6], [1,0.6],\n",
    "              [0,0.4],[0.4,0.4],[0.6,0.4],[1,0.4],\n",
    "              [0,0], [0.4,0], [0.6,0], [1,0]]\n",
    "\n",
    "    splines = [[0,1],[1,2],[2,3], \n",
    "               [0,4],[1,5],[2,6],[3,7],\n",
    "               [4,5],[5,6],[6,7],\n",
    "               [4,8],[5,9],[6,10],[7,11],\n",
    "               [8,9],[9,10],[10,11],\n",
    "               [8,12],[9,13],[10,14],[11,15],\n",
    "               [12,13],[13,14],[14,15]\n",
    "              ]\n",
    "\n",
    "    longs = [0,2,3,4,5,6,7,9,14,16,17,18,19,20,21,23]\n",
    "\n",
    "    surfaces = [[0,4,7,3],[1,5,8,4],[2,6,9,5],[7,11,14,10],[9,13,16,12],[14,18,21,17],[15,19,22,18],\n",
    "             [16,20,23,19]]\n",
    "\n",
    "    for x, y in points:\n",
    "        g.point([x, y])\n",
    "\n",
    "    for i, s in enumerate(splines):\n",
    "        if i in longs:\n",
    "            g.spline(s, el_on_curve=2*el_on_curve)\n",
    "        else:\n",
    "            g.spline(s, el_on_curve=el_on_curve)\n",
    "\n",
    "    for f in surfaces:\n",
    "        g.struct_surf(f)\n",
    "\n",
    "\n",
    "    bcs = {100:[21,22,23],\n",
    "           200:[0,1,2],\n",
    "           300:[3,10,17],\n",
    "           101:[6,13,20],\n",
    "           201:[15],\n",
    "           301:[12],\n",
    "           401:[8],\n",
    "           501:[11],\n",
    "           }\n",
    "\n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "    return g\n",
    "\n",
    "\n",
    "def dir100(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir200(x,y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir300(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def neu101(x,y):\n",
    "    return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def neu301(x,y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def neu401(x,y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def neu501(x,y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "from fem import Poisson2D\n",
    "geometry = create_geometry(el_on_curve=5)\n",
    "Dirichlet = [100,200,300]\n",
    "Neumann = [101,201,301,401,501]\n",
    "hole = [[201,301,401,501]]\n",
    "\n",
    "BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "          101:neu101, 201:neu201, 301:neu301, 401:neu401,501:neu501}\n",
    "\n",
    "\n",
    "model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc)\n",
    "model.mesh(showMesh=True)\n",
    "model.efem()\n",
    "# model.fem()\n",
    "\n",
    "# model.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "f9496677-8681-4b1b-b906-3d3eebc45820",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def create_geometry(el_on_curve=20):\n",
    "    # externel boundary ID\n",
    "    dir_id_1 = 100\n",
    "    dir_id_2 = 200\n",
    "    dir_id_3 = 300\n",
    "    neu_id_1 = 400\n",
    "    # hole boundaryd ID\n",
    "    hole_neu_id_1 = 500\n",
    "    hole_neu_id_2 = 600\n",
    "    hole_neu_id_3 = 700\n",
    "    hole_neu_id_4 = 800\n",
    "\n",
    "    # Geometry definition\n",
    "    g = cfg.Geometry()\n",
    "    g.point([0.0, 0.0], ID=0)\n",
    "    g.point([1.0, 0.0], ID=1)\n",
    "    g.point([1.0, 1.0], ID=2)\n",
    "    g.point([0.0, 1.0], ID=3)\n",
    "    g.point([0.4, 0.4], ID=4)\n",
    "    g.point([0.6, 0.4], ID=5)\n",
    "    g.point([0.6, 0.6], ID=6)\n",
    "    g.point([0.4, 0.6], ID=7)\n",
    "    g.spline([0, 1], el_on_curve=el_on_curve, marker=100, ID=0)\n",
    "    g.spline([1, 2], el_on_curve=el_on_curve, marker=101, ID=1) \n",
    "    g.spline([2, 3], el_on_curve=el_on_curve, marker=200, ID=2) \n",
    "    g.spline([3, 0], el_on_curve=el_on_curve, marker=300, ID=3)\n",
    "    g.spline([4, 5], el_on_curve=el_on_curve/4, marker=201, ID=4)\n",
    "    g.spline([5, 6], el_on_curve=el_on_curve/4, marker=301, ID=5)\n",
    "    g.spline([6, 7], el_on_curve=el_on_curve/4, marker=401, ID=6)\n",
    "    g.spline([7, 4], el_on_curve=el_on_curve/4, marker=501, ID=7)\n",
    "    g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "\n",
    "\n",
    "    bcs = {100:[0],\n",
    "           200:[2],\n",
    "           300:[3],\n",
    "           101:[1],\n",
    "           201:[4],\n",
    "           301:[5],\n",
    "           401:[6],\n",
    "           501:[7],\n",
    "           }\n",
    "\n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "    return g\n",
    "\n",
    "\n",
    "def dir100(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir200(x,y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir300(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def neu101(x,y):\n",
    "    return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def neu301(x,y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def neu401(x,y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def neu501(x,y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def main():\n",
    "    geometry = create_geometry(el_on_curve=20)\n",
    "    Dirichlet = [100,200,300]\n",
    "    Neumann = [101,201,301,401,501]\n",
    "    hole = [[201,301,401,501]]\n",
    "\n",
    "    BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "              101:neu101, 201:neu201, 301:neu301, 401:neu401,501:neu501}\n",
    "\n",
    "\n",
    "    model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc)\n",
    "    model.mesh(showMesh=True)\n",
    "    model.efem()\n",
    "    # print(model.FU)\n",
    "    # model.fem()\n",
    "\n",
    "    # model.display()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a7c985e-022c-4eed-9ed3-e9c0f1e4ed0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mayavi import mlab\n",
    "\n",
    "x = model.Node[:,0]\n",
    "y = model.Node[:,1]\n",
    "z = np.zeros_like(x)\n",
    "triangles1 = model.Element\n",
    "triangles2 = model.Element[list(model.crossElement.keys()),:]\n",
    "\n",
    "mlab.triangular_mesh(x, y, z, triangles1, color=(1,1,1),representation='wireframe' )\n",
    "mlab.triangular_mesh(x, y, z, triangles2, color=(1,0,0),representation='surface' )\n",
    "\n",
    "for cluster in model.NeuCluster:\n",
    "    i,_ = cluster[0]\n",
    "    x = model.Node[i][0]\n",
    "    y = model.Node[i][1]\n",
    "    z = 0    \n",
    "    mlab.points3d(x,y,z,1, color=(0,1,0),scale_factor=0.01)\n",
    "    for i, j in cluster:\n",
    "        x = model.Node[i][0]\n",
    "        y = model.Node[i][1]\n",
    "        z = 0\n",
    "        u = model.Node[j][0] - x\n",
    "        v = model.Node[j][1] - y\n",
    "        w = 0\n",
    "        mlab.quiver3d(x, y, z, u, v, w, scale_factor=1)\n",
    "\n",
    "for edge in model.crossEdge:\n",
    "    i,j = edge\n",
    "    x = [model.Node[i][0], model.Node[j][0]]\n",
    "    y = [model.Node[i][1], model.Node[j][1]]\n",
    "    z = [0,0]\n",
    "    mlab.plot3d(x, y, z, color=(0,0,1),tube_radius=0.005)\n",
    "        \n",
    "        \n",
    "mlab.view(0,0)\n",
    "# print(mlab.view())\n",
    "mlab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6c531ef-1065-4949-bfe2-b55dff73868d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
